program main

 implicit none

 call star ()

 stop
end

!************************************************************

subroutine star ()

 implicit none

 integer (kind=4),parameter :: n=2
 integer (kind=4), parameter :: l=256
 integer (kind=4) k
 external star_f
 real (kind=8) x0,x1,xmax,dx
 real (kind=8) u0(n),u1(n)
 real (kind=8) the0, the1
 real (kind=8) Rg,cp,cap,Rs
 real (kind=8) t00,rh00,M00,PI,g0,g1,Gc

 PI = 3.14159
 Gc = 6.674e-11
 
 !Rs = 1.36206e9          !1.1Myr
 !rh00 = 701.952          !1.1Myr
 !t00 = 2.402835e6        !1.1Myr 

 !Rs = 7.87093e8          !6Myr 
 !rh00 = 3570.057         !6Myr
 !t00 = 4.122538e6        !6Myr 

 !Rs = 6.90278e8          !9.25Myr 
 !rh00 = 5290.174         !9.25Myr
 !t00 = 4.677848e6        !9.25Myr 

 Rs = 6.07751e8          !14Myr 
 rh00 = 8006.324         !14Myr
 t00 = 5.201200e6        !14Myr 

 x0 = 0.1*Rs
 xmax = 0.95*Rs
 dx = (xmax-x0)/(l-1)
 Rg = 13732.
 cp = 2.5*Rg
 cap = Rg/cp

 M00 = 4.*PI*rh00*x0**2  
 u0(1) = t00
 u0(2) = rh00

 the0 = u0(1) * (t00*rh00 / (u0(1)*u0(2)) )**cap
 do k=1,l
  write( *, '(2x,g14.6,2x,g14.6,2x,g14.6,2x,g16.8)') x0,u0(1),u0(2),the0 
  x1 = x0 + dx
  call rk4vec(x0,n,u0,dx,star_f,u1)
  the1 = u1(1)*( (t00*rh00) / (u1(1)*u1(2)) )**cap
 
  x0 = x1
  u0(1:n) = u1(1:n)
  the0 = the1

 end do
 
 return
end

!*************************************************************

subroutine star_f (x,n,u,uprime)

 implicit none
 
 integer (kind=4) n
 real (kind=8) x
 real (kind=8) u(n)
 real (kind=8) uprime(n)
 real (kind=8) g0,g
 real (kind=8) Rg
 real (kind=8) mr
 real (kind=8) m1,m2
 real (kind=8) ms,me, me2
 real (kind=8) width
 real (kind=8) xb,Rs
 

 !Rs = 1.36206e9    !1.1Myr
 !g = 26.           !1.1Myr   Escolha para obter o contraste de densidade proximo ao modelo de evolucao estelar
 !g0 = g*(-0.4745+1.245e-8*x-1.234e-17*x**2+3.177e-27*x**3)
 
 !Rs = 7.87093e8    !6Myr 
 !!g = 80.95         !6Myr
 !g = 73.83         !6Myr FIT
 !g0 = g*(-0.4716+2.147e-8*x-3.649e-17*x**2+1.605e-26*x**3)
 

 !Rs = 6.90278e8    !9.25Myr 
 !g= 104.61        !9.25Myr
 !g= 96.23          !9.25Myr FIT
 !g0 = g*(-0.4565+2.422e-8*x-4.683e-17*x**2+2.339e-26*x**3)
 !xb= 0.20*Rs

 Rs= 6.07751e8    !14Myr 
 !g= 140.07       !14Myr
 !g= 131.19       !14Myr FIT
 g= 132.5         !14Myr Escolha para obter o contraste de densidade proximo ao modelo de evolucao estelar
 g0 = g*(-0.3598+2.571e-8*x-5.692e-17*x**2+3.260e-26*x**3)
 xb=0.31*Rs  
 
 ms= 1.5     !Modelo Isentropico
 me= 1.4998   !Ambiente Convectivo
 me2=1.65      !Ambiente Radiativo
 width= 0.01*Rs 
 
 !mr= ms                                          !Modelo Isentropico
 mr= me2 - (me2-me)*0.5*(1.+erf((x-xb)/width))    !Estado Ambiente para modelo com duas zonas
 !mr= me                                          !Estado Ambiente estrela completamente convectiva

 Rg = 13732.
 
 uprime(1) = -g0 / ((mr+1.)*Rg)
 uprime(2) = (u(2) / u(1)) * (-g0/Rg - uprime(1))
 
 return 
end

!**************************************************************

subroutine rk4vec (t0,m,u0,dt,f,u)

 implicit none

 integer (kind=4) m
 real(kind=8) dt
 external f
 real (kind=8) f0(m),f1(m),f2(m),f3(m)
 real (kind=8) t0,t1,t2,t3
 real (kind=8) u(m),u0(m),u1(m),u2(m),u3(m)

 call f (t0,m,u0,f0)
 
 t1=t0+dt/2.0D+00
 u1(1:m)=u0(1:m)+dt*f0(1:m)/2.0D+00
 call f (t1,m,u1,f1)

 t2=t0+dt/2.0D+00
 u2(1:m)=u0(1:m)+dt*f1(1:m)/2.0D+00
 call f (t2,m,u2,f2)

 t3=t0+dt
 u3(1:m)=u0(1:m)+dt*f2(1:m)
 call f (t3,m,u3,f3)

 u(1:m)=u0(1:m)+dt*(f0(1:m)+2.0D+00*f1(1:m)+2.0D+00*f2(1:m)+f3(1:m))/6.0D+00
 
 return
end

